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Abstract 

A MATLAB program has been developed for simulation of polymer blend 
self-assembly with nano-scaled features. The Cahn-Hilliard equation is imple- 
mented to calculate the free energy profile of the polymer blends. The Flory- 
Huggins type of energy is used to estimate the local free energy. The program 
is capable of quantitatively simulate the phase separation of polymer blends. 
The effects such as the substrate functionalization, solvent evaporation, and 
polymer materials properties are included in the program. The program can 
estimate the model parameters from the real experimental processing param- 
eters and the material properties. The simulation results can be evaluated 
quantitatively and compared with the experimental results with analysis tools 
included in the program. 



Introduction 

decades due to its wide applications |1|] . The numerical simulation of this pro- 
cess can be used to research the mechanisms of phase separation of polymer 
blends and predict the unobservable process states and unmeasurable mate- 
rial properties. The mathematical principles and numerical simulation of self- 



assembly via phase separation has been extensively studied 2 



111]. But few 



specific software toolkit have been developed to efficiently investigate this phe- 



nomenon. 



A computer program is developed in MATLAB for the numerical simulation 
of the polymer blends phase separation. With this software, the mechanisms 
of the phase separation are investigated. Also the mobility, gradient energy 
coefficient energy, and the surface energy in the experiment are estimated with 
the numerical model. The software can evaluate the physical parameters in 
the numerical model by implementing the real experimental parameters and 



materials properties. The numerical simulation results can be analyzed with 
the software and the results from the simulation software can be validated with 
the experimental results. 
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Figure 1: Screenshot of the simulation program graphical user interface. 



Fundamentals 



The numerical model for phase separation of polymer blends are established and 



validated with experimental results work 



12l |. The free energy profile during 



the phase separation in a inhomogeneous mixture is described by Cahn-Hilliard 



Equation 



13l4l6l]. as shown below, 



i^(Ci,C2,C3)= [ 

Jv 



fiCl,C2,C3)+ Y, [Ki{VQf]\dV 



(1) 



1=1,2,3 



where / is the local free energy density of homogeneous material, (pi is the 
lattice volume fraction of component i, and kj is the gradient energy coefficient 
for the component i. The total free energy of the system is composed by two 
items as shown in Equation [TJ The first item is the local free energy and the 
second is the composition gradient contribution to the free energy. 

In our study, the local free energy is in the form of Flory-Huggins equation, 
which is well know and studied for polymer blends 17|] The ternary Flory- 
Huggins Equation is shown as follows. 



/(Ci,C2,C3 



J^ (^ In Ci + -^ In C2 + C3 In C3 

Vsite \mi 1712 



X12C1C2 + X13C1CS + X23C2C3 



(2) 



where R is the ideal gas constant, T is the absolute temperature, VsUe is the 
lattice site volume in the Flory-Huggins model, rrii is the degree of polymeriza- 
tion of component i, and Cj is the composition for the component i. 

There are some parameters in the numerical model which can not be mea- 
sured directly, such as the gradient energy coefficient and the mobility. These 
parameters have to be estimated from the experimental parameters. The gradi- 
ent energy coefficient, k, determines the influence of the composition gradient 



to the total free energy of the domain. The value of k is difficult to measure 



181] to 



experimentally. Though efforts have been made by Saxena and Caneba p 
estimate the gradient energy coefficient in a ternary polymer system from ex- 
perimental methods, few experimental results are published for our conditions. 
Initially, the value of k can be estimated by the interaction distance between 
molecules 



RTo? 



?,v 



(3) 



site 



where a is the monomer size. A modified equation to calculate k considering 



the effects of the composition is reported by Gennes, et al. 



2C|. 



RTa^ 

where the subscript, z, represents component i. 



(4) 



The mobility is estimated from the diffusivity of the components. The mo- 
bility of the polymer blends with long chains can be estimated by the equation 



as follows 



m, 



Mi 



^i J-^mJ^e'^site 



nii 



RT 



(5) 



where rrii is the degree of polymerization as stated before, D^ is the dif- 
fusivity of the monomer, and A'^e is the effective number of monomers per en- 
tanglement length. Because of the scarce experimental data for N^, a more 
generalized form is employed for our study, 
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(6) 



The time evolution of the composition of component i can be represented 

as Ubd-y], 



dt 
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where the subscripts i and j represent components 1 and 2, and 



(7) 



M,, =(1 - Ci)'M, + C- J]; Mj i = 1, 2; i = 1, 2, 3 

^ij = - I^ [(1 - C'^)Cj] Mi + C.CjMg i = 1, 2; j = 1, 2 



(8) 



where Ci is the average composition of component i. To simphfy the solution 
of Equation[71 kh = kj + ks, and K12 = K21 = ks, where Kj is the gradient energy 
coefficient in Equation [31 

For detailed discussion and practical scientific cases with this software can 



be found in our previous works |22 



m- 



The MATLAB Program for Sim- 
ulation of Polymer Phase Separa- 
tion 

Design Principles 

The program is developed in MATLAB m-code. A graphical user interface 
(GUI) is implemented in the program created with MATLAB GUI editor. 
MATLAB is widely used in scientific computation and has many toolkits and 
commonly used mathematical functionalities. But implementing the software 
in MATLAB the efficiency of development is greatly improved. Also, by devel- 
oping the program in MATLAB, the program is cross platform. 



6 



The software is designed for daily usage of simulation and experiment scien- 
tists. The program is light weighted and programmed with high computation 
efficiency so that it can produce significant science results in a common PC. 
It also extensible to a parallel version or implement code to use the high com- 
putation performance of GPU. The GUI is implemented so that the users can 
conveniently input the experiment parameters. The results as well as the user 
settings can be saved and revisited by the program. Also, for better assistance 
to a real productive environment, the simulation model is carefully designed, 
so that the users provide the real processing and material parameters and the 
program will produce quantitative results comparable to experimental results. 
Analytical tools are also provided with the program for post-processing of the 
results. 

Numerical Methods 

To solve the partial differential equation, the discrete cosine transform spec- 
tral method is employed. The discrete cosine transform (DOT) is applied to 
the right hand side and left hand side of Equation [71 The partial differential 
equation in the ordinary space then transformed into an ordinary differential 
equation in frequency space. When the ODE in frequency space then is solved, 
the results are transformed back to the ordinary space. 

Comparing to conventional finite element method, the spectral method is 
more efficient and accurate. This method enabled the program to solve the 
equation in a reasonable computation time to investigate the changes of the 
phase separation during an real time span long enough to observe the phase 
evolution. The spectral method is only applied to the spatial coordinates since 
the time length of the evolution is not predictable. Actually the real time 
for phase evolution is usually one of the major concerns as the result of the 
simulation. 

The DCT takes a considerable portion of the computation time. Especially 



in a 3-diinensional numerical model, the 3-dimensional DCT function with con- 
ventional approach has a complexity of 0{n'^), which can not be practical for 
real application on a PC. To overcome this computational difficulty, the code 
can either be translated to C code embedded in MATLAB m scripts, or a differ- 
ent mathematical approach can be implemented as well. In this program, the 
DCT is calculated from the fast Fourier transform (FFT) which is optimized 
in MATLAB. 

Quantitative Simulation vi^ith Real Experimental Parameters 

Many of previous numerical simulations in the self-assembly with polymer 
blends phase separation are qualitative other than quantitative. The results 
can only be used to provide non-quantitative suggestions to the experiments. 
While this program implemented a numerical model which quantitatively simu- 
lates the experimental results with the real processing and material parameters. 
Most of inputs in to this program can be directly measured and read from the 
instrument or material labels. For some of the physical parameters such as 
K and the mobility, the program can provide a start value from the calcula- 
tion with the theoretical model. The user may need to validate the value by 
comparing the simulation results to the experimental results. Eventually, a 
more accurate estimation can be find with optimization methods by setting the 
difference between the simulation and experiment results as the cost function. 
Besides the parameters in Cahn-Hilliard equation, other effects such as the 
evaporation, substrate functionalization, and the degree of polymerization are 
also implemented with the real conditions. The final results are saved and 
summarized. The characteristic length of result pattern from simulation and 
its compatibility with the substrate functionalization are calculated. These 
numbers can be used to compare with the experimental results. 



Data Visualization and Results Analysis 

When running the program, the message from the software will be output to the 
working console of MATLAB. The messages will show the current state and real 
time results of the simulation. Also, when the simulation is started, the phase 
pattern will be plotted in a real time plot window. Users can set the frequency of 
real time plot and the scale factor on the domain of the contour plot in the GUI. 
The results of the simulation will be saved to a folder designated by the user. 
The real time plot will be saved to the result folder. The quantitative results 
will be saved as several comma separated values (CSV) text files. The result 
folder can be loaded into the analysis toolkit of the program and the user can 
view the assessment values such as the characteristic length, the compatibility 
parameters, and the composition profile wave in depth direction with convenient 
plotting tools. Usually these results such as the composition profile in each 
direction in the domain are difficult to observe in experiment results. 



Examples 



shown in this paper. The results of numerical simul 



To compare the sim- 



To demonstrate the capability of this program, example simulation cases are 

ation have been validated 
with the experimental results in our previous work 
ulated results with a real experimental system, we directed the morphologies 
of polystyrene (PS) / polyacrylic acid (PAA) blends using chemically hetero- 
geneous patterns. More specifically, alkanethiols with different chemical func- 
tionalities were patterned by electron beam lithography, which were then used 
to direct the assembly of PS/PAA blends during the spin coating from their 



mutual solvent [25|] . The experimental conditions are implemented into the nu- 
merical simulation. The effects such as the substrate functionalization and the 
solvent evaporation are involved in the numerical modeling. The parameters 




-clean 

"ClC 



meter, & ra=2.70359e-01 & [min average max] of C! 
meter, & ra-1 .9S960e-Ql & [uin ^erage max] of C 

meter, & ra=l .57754e-01 & [min average max] of C: 

meter, & ra=l .30933e-01 & [rfin average max] of C. 

meter, &. ra=l .120B5e-01 & [min average max] of C. 

meter, & ra=9.B1146B-02 & [min average max] of C: 

meter. & ra=S.7346Se-02 & [m'n average max] of C 
meter, & ra=7 .S793Se-02 fi [pinn a^/erage max] of ( 
it:0. Iteration 11: Time=0.11 sec. & Lz=1.21e-0G meter, & ra=7.1S363e-02 & [min average max] of (_ 



Figure 2: The simulation is running with the real time plot of the current ternary 
phase morphology. 



difficult to measure are acquired with the optimization methods [22 1. 

Sophisticated techniques are required to investigate the composition profile 



in the depth of the polymer film [26l ]. While the numerical simulation results 
can provide the composition profile in each position of the file, the composition 
profile change in depth direction can be easily accessed. To investigate the 
composition wave allow the direction perpendicular to the film surface, a thick 
film is implemented to the numerical simulation. This kind of film is not only 
difficult to fabricate and characterize in experiments, however in the numerical 
modeling, the user only needs to change the mesh grid domain size. The depth 
profiles with different substrate functionalization are shown in Figure [3l where 
\fs\ denotes the surface energy term from the substrate functionalization. This 
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term will be added to the total free energy on the interface of the polymer film 
and the substrate. The initial thickness of the film is 1 mm and decreases to 8 
fj,m due to the evaporation of the solvent. The thickness results are scaled by 
0.5 to fit in the figures. It can be seen that a higher surface interaction force can 
result in a faster substrate directed phase separation in the film. A stronger 
substrate interface attraction force can direct the phase separation morphol- 
ogy near the substrate surface. While with a lower surface energy, the phase 
separation dynamics in the bulk of the film overcomes the substrate attraction 
force. It can be seen that at 30 seconds, the substrate functionalization has 
little effects on the morphology on the substrate surface. Also, the checker 
board structure can be seen near the substrate surface with a higher surface 
energy [111 ]. 



To investigate the effects of a more complicated pattern, a larger domain 
is simulated. The pattern on the substrate applied on the substrate surface is 
shown in Figure HI The substrate pattern is designed to investigate the effects 
of various shapes and contains components such as squares, circles, and dead 
end lines in different sizes. The initial surface dimensions of the model are 
changed to 12fimxl2^m. The initial thickness of the film is 1mm and shrinks 
during the solvent evaporation. The elements in the modelling is 384x384x16. 
The average composition ratio of PS/PAA is changed to 38/62 to match the 
pattern. The result patterns from the simulation can be seen in Figure [5l 

It can be seen that in a larger domain with complicated substrate patterns, 
the attraction factor has to be increased to obtain a better replication. In 
general, the increase of the attraction factor will increase the refinement of 
the pattern according to the substrate pattern. But since the substrate pat- 
tern has geometrical features in different sizes, the attraction factor has to be 
strong enough to force the intrinsic phase separation with unified characteristic 
length to match the substrate pattern in different sizes. This would be the 
main challenge to the replication of complicated patterns. It has been reported 
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Figure 3: The phase separation in a thick film. 
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Figure 4: The substrate pattern with comphcated features. 



by Ming et. al. |25l ] that the addition of the copolymer can improve the refine- 
ment of the final patterns in experiments. The reason is that the PAA-b-PS 
block copolymer will concentrate in the interface of the PS and PAA domains 
in the phase separation, therefore decreasing the mixing free energy. Funda- 
mentally, the addition of the block copolymer increased the miscibility of the 
two polymers. To simulate these phenomena, the Flory-Huggins interaction 
parameter is decreased from 0.22 to 0.1 to increase the miscibility of PS/PAA 
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Figure 5: The effects of complicated substrate patterns. 
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in the modening. The resuh pattern is also shown in Figure [5l in comparison 
to the cases without the addition of block copolymers. It can be seen that 
the refinement of the phase separated pattern is improved by the addition of 
the block copolymer. The Cs values of the phase separation with complicated 
patten are measured and plotted in Figure [H 

A assessment parameter, Cg, the compatibility parameter is introduced to 
evaluate the replication of the morphology to the substrate pattern, where a 
higher Cg value denotes a better replication of the polymer film morphology 
according to the substrate pattern It can be seen in Figure [6] that the Cg value 
for the system with block copolymer is 7.69E-01, which is higher than the 
system without the block copolymer when attraction forces are the same. The 
decrease of the Flory-Huggins interaction parameter increases the miscibility 
of the polymers, which will decrease the miscibility gap of the polymers, as can 
be seen in Equation [2J The two phase at equilibrium will be less concentrated 
in different types of polymer. This is an issue may need to be concerned when 
the interaction parameter of the two polymers is changed. 
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Figure 6: The effects of complicated substrate patterns. 

Conclusion 

A computer program for simulation of polymer self-assembly with phase sep- 
aration is introduced. The program is developed in MATLAB m code and 
designed to assist the scientists in real working environments. The program is 
able to simulate the experiment results quantitatively with real experimental 
parameters. The unmeasurable physical parameters such as the gradient energy 
coefficient and the mobility can be estimated with the program. The program 
provides a graphical user interface and analytical toolkits. This program can 
help the scientists in research in polymer phase separation mechanisms and dy- 
namics with high efficiency, convenience of usage, quantitative results analysis. 
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and validated reliability. 
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